Typeset with jpsj2.cls <ver.l.2.1> 



Full Paper 



A Numerical Model for Brownian Particles Fluctuating in Incompressible Fluids 

Takuya Iwashita 1 *, Yasuya Nakayama 2 , and Ryoichi Yamamoto 1,3 

1 Department of Chemical Engineering, Kyoto University, Kyoto 615-8510 
2 Department of Chemical Engineering, Kyushu University, Fukuoka 819-0395 
3 CREST, Japan Science and Technology Agency, Saitama 332-0012 



We present a numerical method that consistently implements thermal fluctuations and hy- 
drodynamic interactions to the motion of Brownian particles dispersed in incompressible host 
fluids. In this method, the thermal fluctuations are introduced as random forces acting on the 
Brownian particles. The hydrodynamic interactions are introduced by directly resolving the 
fluid motions with the particle motion as a boundary condition to be satisfied. The validity 
of the method has been examined carefully by comparing the present numerical results with 
the fluctuation-dissipation theorem whose analytical form is known for dispersions of a single 
spherical particle. Simulations are then performed for more complicated systems, such as a 
dispersion composed of many spherical particles and a single polymeric chain in a solvent. 

KEYWORDS: Brownian motion, hydrodynamic interaction, thermal fluctuation, particle dispersions, direct 
numerical simulation, long-time tails, polymer dynamics 



1. Introduction 

The dynamics of solid particles dispersed in host flu- 
ids is very complicated. Although computer simulations 
have been extensively used as a tool to investigate those 
systems, obtaining meaningful results is not yet straight- 
forward even for the simplest case where the particles are 
monodisperse spheres and the host fluid is Newtonian. 
The main difficulty comes from the consistent treatment 
of the so-called hydrodynamic interaction (Hf) between 
particles mediated by fluid motions, which are induced 
by the particle's motion. 

The mathematical expression for HI is greatly simpli- 
fied if the following assumptions are made: 1) the dis- 
persed particles are all spherical; 2) HI acting among par- 
ticles is pair- wise additive; and 3) the motions of the host 
fluid instantaneously follows the motions of the particles 
(the Stokes approximation). Then, the HI is expressed 
as a tensor that is a function of the particle's positions 
and velocities, without explicitly dealing with the fluid 
motions. The Oseen and the Rotne-Prager-Yamakawa 
(RPY) tensors are probably the most well known forms of 
such simplified HI functions: the former neglects the size 
of the particle, while the latter takes into account some 
corrections to the particle size. The Stokesian dynamics 
(SD) method 1 - 1 is a widely used numerical method along 
the lines of solving the tensor equations. It is based on 
the Langevin-type equations for particles implementing 
the RPY tensor and the lubrication correction. The lat- 
ter is necessary when the distance between particles is 
small compared to the particle radius. 

Completely different numerical approaches have been 
developed recently, 2-6 ) where the motions of host fluids 
are explicitly solved along with the motion of the parti- 
cles so that the HI is directly computed without using 
the three assumptions described above. We refer those 
approaches as direct numerical simulation (DNS) ap- 
proaches. An apparent benefit of using DNS approaches 
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for particle dispersions is that the long-time behavior of 
particle motions is reproduced accurately. For example, 
the velocity auto correlation function (VACF) of a fluc- 
tuating particle is expected to show a non-exponential 
power-law relaxation- known as the hydrodynamic long- 
time tail-if the thermal fluctuations are taken into ac- 
count as well. 

Recently, we also proposed an efficient DNS scheme 
for particle dispersions called the Smoothed Profile (SP) 
method, 7 ' 8 ^ where the discontinuous boundaries between 
particles and a fluid are smoothed out by using a contin- 
uous profile function, thereby achieving greater compu- 
tational efficiency. The SP method has been successfully 
applied to some problems, such as the electrophoresis 
of charged colloidal particles; 9 ' however, the effects of 
thermal fluctuations (TF) are neglected. This is not a 
bad approximation for large/heavy particles, but it is 
insufficient for particles whose size is on the order of a 
micrometer or smaller; the coupling of TF and HI be- 
comes crucial in these systems. Owing to some new ex- 
perimental techniques that enable direct examination of 
the properties of Brownian particles fluctuating in a host 
fluid, several interesting phenomena have recently been 
reported, including the non-diffusive behavior of Brow- 
nian particles 10 ' and the rotational-translation coupling 
of a pair of spherical particles 11 ' where the coupling of 
TF and HI plays an essential role. 

In the present study, our main goal is to implement TF 
and HI consistently within a DNS framework for particle 
dispersions. We aim to introduce TF to our original SP 
method to achieve this end. Since theoretical analyses 
have been well established for dilute dispersions, we first 
perform simulations for a dilute dispersion composed of 
a single spherical particle in a cubic box with the peri- 
odic boundary condition. The numerical results for the 
VACF are then compared with analytical solutions based 
on the fluctuation-dissipation theorem. After the valid- 
ity of the method has been confirmed for this simple 
system, simulations arc performed for dense dispersions 
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composed of many spherical particles for which analyt- 
ical solutions are unknown. We furthermore performed 
DNS simulations to examine the dynamics of a single 
polymeric chain fluctuating in a good solvent. The time- 
correlation functions are calculated for each Rouse mode 
of the chain and compared with the predictions of the 
Rouse and Zimm models. 

2. Simulation Method 

We now present our implementation of TF based on 
the Langevin type approach. The equation governing the 
solvent with a density pf and a shear viscosity 77 is a 
modified Navicr-Stokes equation 



Pf{dtv + v ■ Vu) = — Vp - 
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with the incompressible condition V • v = 0, where 
v(x,t),p(x, t) is the velocity and pressure field of the 
solvent. A smooth profile function < 4>(x,t) < 1 dis- 
tinguishes between fluid and particle domains, namely 
<f> = 1 in the particle domain and <f> — in the fluid 
domain. Those domains arc separated by thin intcrfacial 
regions whose thickness is characterized by £. The body 
force 4>f p is introdued to ensures the rigidity of particles 
and the non-slip appropriate boundary condition at the 
fluid/particle interface. The mathematical expressions of 
4> and (j>fp are given in earlier papers 7,8 ' in detail. 

We consider dispersions composed of N p spherical par- 
ticles with a radius a. The motion of the ith particle is 
governed by Newton's equations of motion with stochas- 
tic forces: 
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where Ri, Vi, and fi, are the position, translational ve- 
locity, and rotational velocity of particles, respectively. 
Mi and J, are the mass and the moment of inertia, and 
F^ and N- 1 are the hydrodynamic force and torque ex- 
erted by solvent on the particle. 7,8 ' Ff is direct intcrpar- 
ticle interaction such as Coulombic or the Lcnnard-Jones 
potential. In the present study, we used the truncated 
Lcnnard- Jones interaction: 



U L j(n 3 ) = 
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\Ri — Rj\. The parameter e characterizes 
the strength of interactions, and a = 2a represents the 
diameter of particles. Gf and Gf are random forces 
and torques due to thermal fluctuations, which has the 
following properties 



(G?(t))=0, (G?(*)G7(0)> = a n I5(t)Si 



(5) 



where the square brackets denote taking an average over 
an equilibrium ensemble. a n represents the noise inten- 
sity for each degree of freedom of the translation (n = V) 
and rotation (n = fi) of the particles. Each noise in- 
tensity is controlled so that the variance of the transla- 
tional and rotational velocity has a constant value; that 



(Vi' 



C\ and (ttj) = C2, where C\ and C2 are 



constant numbers. The time evolution of the noise inten- 
sity is described by a v (t + At) = a v (t)e 1 - v - 2 ( t '/ Cl and 
a n (t + At) = a !! (t)e 1_fi i (t)/C2 , where At is the discrete 
time interval which plays a role of thermostat. 

The temperature of the system is defined by the diffu- 
sive motion of the dispersed particles. The translational 
particle temperature ksT v is determined by the long- 
time diffusion coefficient D v of a spherical particle in the 
Stokes-Einstein relation ksT v = QnrjaD v where D v is 
obtained from computer simulation. Similarly, the rota- 
tional particle temperature ksT n can be determined by 
the rotational diffusion coefficient D n . 

There are several advantages to the Langevin approach 
compared with the fluctuating hydrodynamics approach, 
for which TF is introduced in the Navier-Stokes equa- 
tion as stochastic stresses to be defined on N d grid 
points of fluid simulations. One important advantage is 
the computational efficiency: while a -/V d spatial grid re- 
quires generating 0(N d ) random numbers for the fluctu- 
ating hydrodynamics, our method requires 0(N p ) ran- 
dom numbers for a dispersion composed of -/V p (<C N d ) 
particles. Second, if we consider the solvent as a complex 
fluid-with arbitrary constitutive equation-, our method 
has another merit: in the fluctuating hydrodynamics ap- 
proach, the friction tensor for complex fluids is required, 
but not here. Derivations of the friction tensor for com- 
plex fluids often have theoretical difficulties, and even 
if the friction tensor is obtained, there may be a large 
computational cost. 

3. Test of the Simulation Method 

In order to test our simulation method, we consider 
the translational motion of a particle dragged with a 
constant external force Fq in a Newtonian fluid. With 
the force turned on, a dragged particle and the solvent 
have a steady state solution: the particle has a constant 
velocity along a;-direction until t = 0. Then, at t > 0, the 
external force is removed and the particle and solvent 
relax toward a rest state due to dissipation. 

Numerical simulation has been performed in three di- 
mensions with periodic boundary conditions. The lat- 
tice spacing A is taken to be the unit of length. Other 
units are defined so that we can set r\ = 1 and pf = 1 
in cq.(l), namely the units of time, mass, pressure are 
given by pfA 2 /i], p/A 3 , and rj 2 /pfA 2 , respectively. For 
a spherical particle, the radius a = 5, the thickness £ = 2 
and the particle density p p = I was used. In Fig. I, we 
plotted —ii(t)/F for system size 32 3 , 64 3 and 128 3 where 
u(t) denotes the velocity in the x-direction of a dragged 
particle. The response function almost coincides with the 
analytical solution based on a time-dependent friction 12 ' 
when the system size is larger. The analytical solution 
is provided in detail in Appendix. For the long-time re- 
gion, the regression of the velocity shows a power-law 
decay Bt~ 3 ^ 2 where B — I/12p/(7ri//) 3 / 2 and the kine- 
matic viscosity Vf — T]/pf, which depends only on the 
hydrodynamic property of the solvent. 

A basic relation between the relaxation response of a 
dragged particle and the velocity autocorrelation func- 
tion of a Brownian particle is the fluctuation-dissipation 
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Fig. 1. The response functions, the translational velocity of a 
dragged particle at a constant external force Fq, for system size 
32 3 , 64 3 and 128 3 . The response function is affected by the finite- 
size effect. Dotted line is the analytical solution of a drag prob- 
lem with the time-dependent friction. 12 ' Dash-dotted line is the 
algebraic power Bt~ 3 / 2 where B = l/12pf (irvf) 3 / 2 . 
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~Fq dt 



= ^{Vi{t).VM) 



(6) 



where (Vi(t) ■ Vi(0))/3 is the translational velocity auto- 
correlation function (VACF) for a Brownian particle and 
Pv = l/k B T v . The relation holds that the response 
function of a dragged particle with external force Fq is 
equal to the VACF of a particle in thermal equilibrium. 
Under the same computational conditions as the relax- 
ation experiment, the VACF in thermal equilibrium has 
been calculated at the volume fraction $ = 0.002 and 
system size 64 3 . Figure 2 shows the VACF for a Brownian 
particle. We can calculate the long-time diffusion coeffi- 
cient D v by the mean-square displacement of a Brownian 
particle; that is, lim t _« X) (|.Rj(*) - -R 4 (0)| 2 ) = 6D v t. 
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Fig. 2. Pluses(+): Translational velocity autocorrelation function 
for Brownian particles at fiy = 1-2 and $ = 0.002. The solid line 
is the analytical solution of a drag problem for the translational 
particle velocity. The dotted line is the algebraic power £?4~ 3 / 2 
where B = l/12pf(iTUf) 3 / 2 . The dashed line is the exponential 
decay obtained from Markovian Langevin equation. 



When simulating a Brownian particle with the HI, 
the diffusion coefficient is affected by finite-size effects 
and is given by D v = k B T v /6irr)aK($) where K(<Z>) 



represents the effect of the periodic boundary condi- 
tion. 13,14 ' Taking the finite-size correction into account, 
the diffusion coefficient at infinite dilution is obtained 
as Djnf = D v K($). The translational particle temper- 
ature is estimated to be k B T v = 6irr)aDY n f — 0.83. In 
Fig. 2, (iy (Vi(t) • Vi(0))/3 approaches to the analytical so- 
lution in the long-time region, shows the power-law decay 
Bi~ 3 / 2 , and gives the response function of a dragged par- 
ticle. In the short-time region, a gap between simulation 
and the analytical solution is observed. The amplitude 
of a particle's velocity is related to the equipartition law 
of energy; that is, (V?) = 3k B T v /M e where M e is the 
effective mass. The effective mass can be obtained via a 
hydrodynamic calculation as M e = Mi + 0.5mo, where 
mo = Airpfd 3 /3 is the mass of fluid displaced by a spheri- 
cal particle with radius a. In this simulation, the effective 
mass is estimated to be M e ~ 4.2M,, which is notably 
greater than the hydrodynamic effective mass. The dis- 
agreement may be due to the influence of the artificial 
smoothed boundary used between a particle and fluid. 
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Fig. 3. The spring constant dependence of the 
limt— >oo{Ar 2 (t))fc/2. kgT v is the translational particle 
temperature, which is obtained by the diffusion coefficient for a 
Brownian particle. 



Furthermore, we investigate Brownian particles in har- 
monic potentials. The harmonic potential with a spring 
constant k is introduced by adding Ff x = —k(Ri~R^ q ), 
where is its equilibrium position, to the equations 
of motion of particles. The mean-square displacement of 
a Brownian particle trapped in a harmonic potential is 
given by: 

2k B T 



lim (Ar 2 {t)) = 



(7) 



k ' 

where (Ar 2 (t)) = (\Ri(t) - il l (0)| 2 )/3. In the simulation 
for the particle number N p = 8, each particle is trapped 
at the grids of the fee lattice by harmonic potentials and 
the direct interaction between particles is ignored. For 
k = 0.5, 1.0,2.0 and 4.0, the limt^ oc (Ar 2 (t))fc/2 is plot- 
ted in Fig. 3. The result approaches k B T v for k — > and 
is consistent with the results obtained from the diffusion 
coefficient for a Brownian particle. 

We investigate the response function of the rotational 
motion of a particle with a constant external torque 
N$ in a similar fashion as the investigation of the par- 
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tide's translational motion. Figure 4 shows the relax- 
ation response —uj/Nq where cu(t) denotes the particle's 
rotational velocity. The response function almost coin- 
cides with the analytical solution, and a power-law de- 
cay Ct -5 / 2 is observed in the long-time region where C = 
71-/32,0/ (tw/) 5 / 2 . The FDT for the rotational velocity of a 
particle is also investigated. The rotational velocity au- 
tocorrelation function (ni(t)-{li(0)}/3 (RVACF) is given 
in Fig. 4, where fij(t) denotes the rotational velocity of 
a Brownian particle. In the Grccn-Kubo formula for the 
rotational velocity of the particle, the rotational particle 
temperature is estimated to be ksT = 8iTj]a 3 D n ~ 1. 
The RVACF also shows a power-law decay Ct~ 5 ^ 2 in the 
long-time region. The effective moment of inertia I e can 
be estimated by the value of the same time correlation 
as J„ ~ 2.6J,. 
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Fig. 4. Crosses(x): Rotational velocity autocorrelation function 
for Brownian particles at /3q = 1 and <t> = 0.002. Pluses(+): The 
relaxation response of rotational velocity of a particle with a con- 
stant external torque Nq. The solid line is the analytical solution 
of a drag problem for the rotational particle velocity. The dashed 
line is the algebraic power Ct~ 5//2 where C = 7r/32p/(irv/) 5 ' /2 . 
The dashed-dotted line is the exponential decay obtained from 
Markovian Langevin equation. 



4. Applications 

As a demonstration of a DNS incorporating thermal 
fluctuation of particles, our method is applied to a many- 
particle system and a dilute polymeric chain. 

4-1 Many particles system 

In Fig. 5, the translational velocity autocorrelation 
function for each volume fraction $ is presented. As the 
volume fraction is increased, it is found that the relax- 
ation was more rapid than in low volume fractions, and 
the power-law long-time tail of the VACF gradually dis- 
appear in Fig. 5(a). Figure 5(b) shows that for $ > 0.4, 
the VACF has a negative overshoot, which represents an 
oscillative motion of a tagged particle due to a transient 
cage composed of surrounding particles. Compared with 
Brownian dynamics without HI, the decay of the corre- 
lation for a high volume fraction is much slower with HI, 
probably due to the lubrication force between particles. 
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Fig. 5. The translational velocity autocorrelation function for 
each volume fraction. The parameters of the simulation are 
pf = l,n = l,p p = 1, o = 5,§ = 2, j3y = 1 and /3q = 1. 
Fig. 5(a) is the log-log plot. Fig. 5(b) is the normal plot. 



4-2 A dilute polymeric chain in good solvent 

The role of the HI is important in the dynamics of a 
dilute polymeric chain in a good solvent, and the simple 
theoretical model is known as the Zimm model. 15 -' This 
model treats the HI between beads as a hydrodynamic 
mobility matrix, such as the Oseen tensor. Some groups 
have studied a single polymeric chain with the HI using 
similar hybrid simulation methods, and their results are 
in agreement with Zimm theory. 3 ' 16 ' 17%> Here, we reex- 
amine the validity of the Zimm model using our present 
DNS method since it is supposed to be more accurate 
than other methods used previously. 

As a model of a polymeric chain, we study a bead- 
spring model with a truncated Lennard-Jones potential 
and a finitely extensible non-linear elastic (FENE) po- 
tential 18 ' : 

U F (r) = -h c R 2 Q Hl - (r/i? ) 2 }, (8) 

where k c = 30e/a 2 , R$ = 1.5cr and r is the distance 
between the neighboring beads. The position vector of 
a bead is described by R n {t) with n = 0, 1, N c h — 1 
where N c h denotes the total number of beads. 

The static property of a polymeric chain is charac- 
terized by the static exponent which is defined as 
(Rq) oc N^b 2 for large N c h, where Rq is the radius 
of gyration and b is the average bond length. The static 
exponent is related to the size of a polymeric chain, which 
is v = 0.5 for a Gaussian chain and v ~ 0.6 for a self- 
avoided chain. The static exponent v of a polymeric chain 
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can be calculated via the static structure factor S(k). 
Figure 6 displays the static structure factor for the bead 
numbers N c h = 10, 15. In the range of Rq 1 <C k -C b^ 1 , 
S(k) obeys the scaling relation S(k) cx fc -1 / 1 ' and can de- 
termine the static exponent v ~ 0.62 by fitting a power 
law to our data. 




Fig. 6. The static structure factor of a chain, for N c ^ = 10 
(Pluses), 15 (Crosses). The parameters of the simulation are 
Pf = 1,77 = X,pp = 1, a = 5, £ = 2, f3y = 1 and (3q = 1. 
System size is 64 3 grids. 



To analyze the relaxation dynamics of a polymeric 
chain, the real space motion is decomposed into a set 
of the Rouse modes (p — 0, 1, • • • , N c h — 1), 



X p (t) = — V R n (t) cos[(n - l/2)pn/N ch }. 

J * ch 



(9) 



Within the approximations of the Zimm theory, the au- 
tocorrelation function of the Rouse mode decays expo- 
nentially as (X p (t) ■ X p (0))/(Xp) = cxp(— t/rp), where 
t p denotes the relaxation time of the Rouse mode. The 
Zimm theory predicts the relation between the static 
exponent v and the relaxation time r„. The predic- 



tion is t 

r.2 



for the continuous model, and t„ 



p - P 

p~ sin -2 (pn /2N c h) for the discrete model. 

Figure 7 shows the normalized autocorrelation func- 
tion of the Rouse mode for N c h = 10. The relaxation 
times are obtained by a fitting in the exponential short- 
time regime t € [50 : 1000]. Figure 8 shows the mode (p-) 
dependence of the relaxation time. By fitting for p < 5, 
the p-dependence of the relaxation time is estimated as 
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Fig. 7. Normalized autocorrelation functions of the Rouse mode 
X p for various mode p. The chain length is N c h = 10. 
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Fig. 8. Rouse mode p-dependence of the relaxation time for 
chain length N c h = 10. The dashed line is a fitting power 
law to simulation data for p < 5. The theoretical prediction 
of the discrete Zimm model is indicated as a dotted line, t p ~ 



paring the present numerical results with the fluctuation- 
dissipation theorem for a dispersion of a single spheri- 
cal particle. Simulations are then performed for disper- 
sions of many spherical particles, and also for a polymeric 
chain in a fluid. In the former case, we found that the hy- 
drodynamic long-time tail in the VACF-clcarly observed 
for a single particle dispersion-becomes weak with in- 
creasing volume fraction of the particles. In the latter 
case, we found that our numerical results coincide quite 
well with the theoretical predictions of the Zimm model. 
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In Fig. 8, the prediction of the discrete Zimm Appendix: A Dragged Particle with the Time- 
dependent Friction 

From the hydrodynamics, we can obtain the time- 
dependent friction 12 ' 



model for the p-dependence of r p is also plotted using v 
obtained from S(k). The numerical results show a good 
agreement with the prediction of the Zimm model. 



5. Summary 

We have developed a numerical method for consis- 
tently implementing thermal fluctuations and hydrody- 
namic interactions into models of the motions of Brown- 
ian particles dispersed in incompressible host fluids. We 
represented the thermal fluctuations by random forces 
acting on Brownian particles and the hydrodynamic in- 
teractions by directly resolving the fluid motions. The va- 
lidity of the method has been examined carefully by com- 
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The first term is the standard Stokes resistance. The sec- 
ond term represents the additional mass, which is related 
to the acceleration of the particle. The third term rep- 
resents the memory effect, which addresses the temporal 
decay of the fluid's momentum. 

The equation of motion of a dragged particle with a 
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time-dependent frictional force £ and a constant external 
force Fq is described by Mu{t) = ((t) + Fq, where M is 
the mass of a particle. After the particle has reached 
a steady state (a constant velocity), the external force 
is removed. Then, the regression of the velocity can be 
written as 



du(t) 
dt 



F 



dy voy/ye 



M eff J Q 7T |l-y|2 +(To V (A ' 2) 

where er = (%/(2p p + p/)) 1/2 , M eff = M + 0.5m , 
p p the particle density and tb = M/6Trr]a. From the 
fluctuation-dissipation relation 3k ^ T v {Vi[t) ■ Vi(0)) = 
— j^du(t)/dt, a hydrodynamic velocity autocorrelation 
function of a Brownian particle can be calculated exactly. 

Similarly, the regression of the rotational velocity lo of 
a spherical particle with a moment of inertia / can be 
written as 



duj(t) 
dt ' 

N 



1 f°° dy 

— / ^- exp(-yt/r / ) 
T f J Sir 



„3/2 



8nr]a 3 Tf 

where No is a constant external torque, 77 = a 2 /vj and 
7> = I/8nr)a 3 
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